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INTRODUCTION 


This  report  present!  the  result!  of  in  analytical  study  for  predicting  nuclear  blast  wave  propagation 
in  air  entrainment  systems  of  hardened  facilities.  This  is  the  final  report  of  the  study,  and  the  information  that 
is  presented  supplements  the  initial  report  titled  Shock  Wave  Propagation  Through  Air  Fntrainmcnt  Systems— 
Phase  I | ' 1 •*  The  analytical  approach  to  the  problem  is  discussed  in  the  initial  report  with  pre’iminarv  results 
frem  two  computer  codes.  These  initial  computer  codes  differed  basically  in  .he  finite-difference  method  that 
was  utilized  for  solution  of  the  fundamental  gas  dynamics  partial  differential  equations.  One  rode  used  the 
pseudo-viscosity  method  |2|  in  a Lagrange  formulation,  and  the  other  used  a modified  Lax-Wendroff  two-step 
(2|  from  S.  Burs.ein  1 3. 4|  in  an  Kulcrian  formulation.  The  solutions  presented  in  the  inir;*!  report  are  for 
oae-dimensional  shock  wave  propagation  in  a constant  area  duct  with  the  effects  of  visccsity  at  the  duct  wall  - 
included. 

Since  the  initial  report  was  issued,  the  computer  codes  have  been  modified  considerably.  The  l.agr-.nge 
computer  code  has  been  altered  to  include  a variable  cross-scctiona*  area  capability  and  wave  propagation  through 
a junction  of  three  ducts.  The  Kulcrian  computer  code  has  been  changed  to  the  original  Lax-Wcndroff  two-step 
scheme.  A modification  has  also  been  made  to  the  Kulcrian  code  to  include  pscuc'.o-viscosity  in  the  finite- 
difference  equations  for  improved  stability  characteristics  as  recommended  in  References  2 and  4.  Both 
computer  codes  have  been  refined. 

Following  is  a description  of  the  I jy range  and  Eulerian  computer  codes,  including  the  background 
for  the  finite-difference  equation  theory  and  input  and  output  formats.  The  version  of  the  Kulcrian  code 
is  for  shock  propagation  with  maximum  air  temperatures  of  1,000°K.  The  Lagrange  code,  however,  includes 
equations  of  state  that  are  appropriate  for  air  temperatures  up  to  24,000°K  and  can  be  used  for  nuclear  blast 
wave  calculations. 


LAGRANGE  COMPUTER  CODE  DESCRIPTION 

The  Lagrange  computer  code  is  a one-dimensional  variable  area  code  which  ’ tcludcs  viscous  effects 
(at  a wall).  The  code  is  suitable  for  computation  of  time-dependent  flows,  including  normal  shock  wives,  in 
a single  duct  in  which  the  cross-sectional  area  is  either  constant  or  variable.  Figure  I,  and  in  the  multiple  duct 
configuration  shown  in  Figure  2.  The  basic  equations  and  finite-difference  equations  used  differ  from  those 
given  in  Reference  1 due  to  the  inclusion  of  t.riable  cross-sectional  area.  Two  types  of  boundary  conditions 
can  be  specified  at  the  duct  system  inlet:  a side-on  type  entrance  with  a surface  shock  wave  specified  that 
passes  over  the  entrance  (see  Figure  I ),  or  simulation  of  a shock  wave  originating  at  the  duct  inlet.  In  both 
cases  the  shock  wave  can  be  of  constant  strength  or  a wave  with  exponential  pressure  decay  behind  the  shock 
front.  For  the  configuration  of  Figure  2,  a nuciear  wave  can  be  specified  on  the  surface.  The  boundary  condi- 
tion at  the  exit  of  a duct  is  specified  as  a rigid  wail  except  at  a T-junction. 


Numbers  in  brackets  refer  to  References. 


(a)  Constant  area  duct  with  side-on  entrance. 
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(Ii)  Duct  with  area  expansion  and  side-on  entrance. 

figure  I.  Single  duct  configurations. 
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Figure  2.  Typical  air  entrainment  system. 


Basic  Equations 


The  basic  differential  equations  arc  written  in  the  Lagrange  formulation  with  the  Lagrange  variable 
defined  as  j 


Of 


where  A is  the  duct  cross-sectional  area,  V is  the  volume  per  unit  ’.nass,  and  x is  the  Eulerian  distance  coordinate. 
The  area  A is  considered  a function  of  x.  A relationship  for  V follows  from  Equation  1 that  has  the  form 
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(2) 


The  above  equation  with  the  definition  for  the  velocity  in  Lagrange  coordinates,  u = 3x/3t  yields  the  equation 
for  conservation  of  mass,  including  variable  area,  as 

3V  3u  u V 

aT  " v 37  + X 
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The  inclusion  of  variable  area  doc*  not  alter  the  form  of  either  the  momentum  equation  or  the  energy 
equation.  These  equations  are.  therefore,  the  same  as  for  a constant  area  duct,  a*  gf.cn  in  Kcference  1. 
The  momentum  equation  is 


jnd  the  energy  equation  is 
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where  p is  the  pressure,  e the  internal  energy  per  un>t  mass,  f a wall  friction  coefficient,  and  D the  duct 
diameter.  The  system  of  equations  is  completed  by  rhe  equation  of  state  that,  for  the  analysis  presented  here, 
take*  the  form 


e 


(6) 


In  the  above,  y is  the  adiabatic  exponent  for  a real  gas.  The  determination  of  the  value  of  y is  explained  in 
the  section  on  equations  of  state.  For  purposes  of  computer  computation,  Kquations  2,  4,  5,  and  6 with  the 
definition  for  the  velocity  u are  written  in  fin  tc-difference  form,  utilizing  the  pseudo-viscosity  method  |2| 
of  shock  wave  treatment. 


The  Finite-Difference  Kquations 

The  finite-difference  equations  for  the  variable  area  computer  code  differ  slightly  from  those  for  the 
constant  area  code  reported  in  Reference  1.  The  system  of  equations  for  the  variable  area  case  is  given  below. 
In  these  equations  superscripts  represent  time  increments  and  subscripts  represent  space  increments,  so  that 
u"  denotes  the  velocity  of  the  jlh  grid  point  of  the  finite-difference  mesh  at  the  n,h  time  increment.  Frac- 
tional space  increments  represent  centered  or  mesh  midpoint  values  as  .hown  in  Figu*.  a.  T5v*  finite-difference 
equations  are: 
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In  the  above  q is  the  artificial  viscosity  ami  I*  » p ♦ q.  The  time  step  At"* 1 12  is  the  average  of  previously 
ami  newly  calculated  time  steps,  with  .it1'  * ' being  the  newly  calculated  time  step.  The  symbols  a(  and  a, 
•ire  constants. 

The  stability  criterion  that  must  be  satisfied  for  the  finite-difference  equations  to  be  stable  is.- 
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The  Boundary  Conditions 


Inlet  Boundary  Conditions.  In  the  Lagrange  code  two  types  of  duct  inlet  boundary  conditions  can 
be  specified:  a side-on  type  entrance  or  simulation  of  a shock  wave  originating  at  a duct  inlet.  In  a Lagrange 
formulation  it  is  only  necessary  to  specify  the  pressure  at  the  duct  inlet  and  the  initial  state  of  the  gas  in  the 
duct  to  produce  a shock  wave  in  the  duct.  Values  of  other  parameters,  such  as  velocity,  density,  and  internal 
energy  behind  the  shock  wave,  arc  determined  by  the  conservation  equations  and  equation  of  state.  Since 
the  zones  of  the  finite-difference  grid  arc  not  stationary  in  this  formulation,  a new  zone  needs  to  be  established 
periodically  at  the  ducr  entrance  to  allow  for  mass  inflow.  In  establishing  this  new  zone  two  state  parameters 
have  to  be  specified;  in  this  ease  the  state  parameters  used  are  the  pressure  and  internal  energy.  When  a shock 
wave  of  constant  streng'h  is  established  at  the  duct  entrance,  the  internal  energy  and  pressure  for  a new  zone 
is  determined  uniquely  by  the  shock  Mach  number.  With  the  inlet  shock  Mach  number  specified,  the  initial 
inlet  pressure  is  determined  from  the  relation 

p(0.0)  » pjl  ♦ p’)Ms2  - /rl  W 


where  pffl.O) 


pressure  at  x = 0 for  t = 0. 


Mh  * shock  .Mach  number 

pa  * ambient  pressure 

ft2  = <7  - l>/<7+  I) 


When  specifying  a shock  of  constant  strength,  Kquation  9 is  used  for  t ? 0. 

For  shock  waves  with  pressure  decay  following  the  shock  front  the  duct  inlet  pressure  for  t > 0 is 
aoproximated  by. 


p«>,t) 
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-t/tj 


(10) 


where  t;  is  the  initial  slope  time  intercept  for  the  wave.  This  simple  form  for  the  inlet  pressure  is  adequate 
for  simulating  shock  waves  that  do  not  have  a significant  negative  pressure  phase,  e.g.,  TNT  blast  waves 
outside  of  the  fireball  region  |5|  and  experimental  shock  tube  waves  after  rarefaction  catchup  has  occurred 
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1 1, 6| . To  simulate  nuclear  bfaat  wave*  where  a significant  negative  pressure  phase  occurs  Equation  10  is 
replaced  by  more  appropriate  relationships  such  as  are  given  in  Reference  7.  These  relationships  for  a 
surface  wave  that  are  used  in  side-on  entrance  calculations  are  given  below: 
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n-t,)/D; 

shock  arrival  time 

duration  of  positive  pressure  phase 

shock  peak  overpressure  at  t » ts 


The  quantities  A,,  A2,  Aj,  b,,  b2,  and  bj  are  constants  for  which  values  are  given  in  Reference  7 for  a 1-Mt 
nuclear  burst.  In  addition  to  a relationship  fo»  -he  surface  pressure,  pv  relationships  for  the  dvnamic  pressure, 
Q,,  and  temperature,  are  required. 
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where  w 
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Q,. 


duration  of  positive  velocity  phase 
peak  dynamic  pressure  ar  t » t, 
•hock  temperature  at  t » t% 


Thr  quantities  A4,  Aj,  b4,  bj,  and  b6  are  constants  for  which  values  are  given  in  Reference  7.  Equations  1 1, 
12,  and  1 3 completely  define  the  surface  conditions  in  the  air  above  an  inlet  from  which  other  variables  can 
be  obtained.  The  total  enthalpy  is  required  for  side-on  entrance  simulation  and  is  defined  as 


The  density  p,  is  computed  using  Equations  1 1 and  1 3,  and  the  equation  of  state. 


(14) 
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The  parameter  Z,  is  a function  of  temperature  and  density  and  is  determined  as  explained  in  the  section  on 
equations  of  state.  The  internal  energy,  e,,  is  determined  from 
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where  yt  it  also  a function  of  temperature  ami  density.  The  particle  velocity,  u,,  is  computed  using 
Equation  12  and  the  definition 
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Simulation  of  a side-on  type  entrance  is  achieved  by  using  Equations  9 and  10  or  1 1 to  compute  the 
surface  pressure,  pv  (point  s on  Ki'-urc  2)  and  then  calculating  the  duct  inlet  pressure  (point  e on  Figure  2) 
from  the  empirical  relation, 


Pc 
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(18) 


The  above  relationship  was  obtained  from  data  in  Reference  8.  (In  the  computer  program  Hquation  18  is 
input  in  table  form.)  When  a new  zone  is  formed  at  the  duct  entrance,  the  internal  energy  of  the  new 
zone  is  determined  by  assuming  the  total  enthalpy  at  points  s and  e are  equal.  This  assumption  yields, 


(19) 


where  h(s  is  the  total  enthalpy  at  point  s.  The  equation  of  state  has  been  used  in  obtaining  Equation  19.  The 
value  of  h(1  is  known  since  values  for  all  shock  parameters  are  known  at  point  s;  it  is  given  by  Equation  14. 

The  value  of  ur  is,  however,  unknown  and  is  approximated  by  the  value  calculated  during  the  previous  time 
cycle.  An  iteration  technique  could  be  used  to  improve  the  value  for  ur,  but  a comparison  of  computed  results 
with  experimental  data  showed  this  to  be  unnecessary  (Reference  1 ).  The  pressure  in  the  new  zone  is  initially 
assumed  as  the  mean  between  the  inlet  pressure.  pc,  and  the  pressure  in  the  second  zone.  The  pressure,  internal 
energy,  and  velocity  completely  define  the  initial  state  for  the  new  zone.  This  approximate  method  for  estab- 
lishing a new  inlet  zone  to  allow  mass  inflow  yields  a result  for  the  shock  peak  pressure  that  is  approximately 
5"  • high  and  a shock  speed  that  is  approximately  2%  high.  Rezoning  methods  will  be  discussed  in  more 
detail  in  a separate  section  below. 

Some  discussion  is  in  order  concerning  the  appropriateness  of  Equation  18  in  defining  the  flow  losses 
at  the  duct  inlet.  Equation  18  defines  the  static  pressure  change  through  the  inlet  and  agrees  essentially  with 
data  from  other  sources  than  BRI.,  c.g..  TRW  Systems  investigation.  IITRI  investigation,  and  the  B.E.I..  Decklcr 
and  D.  II.  Male  investigation.  References  9,  10,  and  1 1,  respectively.  Using  the  second  law  of  thermodynamics 
an  expression  characterizing  the  inlet  loss  can  be  derived  if  y and  Z are  assumed  constant  from  point  s to  point 
c (F'gure  2).  This  loss  is  measured  by  the  entropy  change  and  is  given  by 


(20) 


where  S^  and  Sc  arc  the  entropy  at  point  s and  point  e,  respectively,  and  Z and  y are  the  values  at  point  s.  An 
.alternate  form  for  Equation  20  in  terms  of  the  stagnation  pressures  at  points  s and  c is 
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(21) 
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Examination  of  Equation*  20  and  21  clearly  show*  that  the  flow  lot*  through  an  inlet  or  junction  is  determined 
solely  by  the  stagnation  pressure  change;  specification  of  static  pressure  change  only,  as  by  Equation  18,  neglects 
the  density  change  effects.  For  consign  normal  sh-^k  waves  specification  of  the  static  pressure  above  is  correct 
since  all  other  parameters,  such  as  density  and  stagnation  pressure,  are  determined  uniquely  from  the  static 
pressure  behind  the  shock  wave.  However,  for  shock  waves  that  have  a rarefaction  behind  the  shock  surface 
or  in  which  the  air  behind  the  shock  surface  is  heated  by  radiation,  e.g.,  a nuclear  wave,  the  density  change 
should  be  included  as  indicate  J in  Equation  20.  Data  currently  available  do  not  allow  for  inclusion  of  densirv 
effects  except  for  a limited  amount  of  data  in  Reference  9.  To  specify  the  losses  at  an  inlet  more  accurately 
than  is  done  by  Equation  18  additional  data  arc  required  in  which  the  surface  density  (or  a related  state  variable 
such  as  temperature)  is  varied  independently  of  the  surface  pressure. 

T-Junction  Boundary  Condition.  The  T-junction  (Figures  2 and  4)  for  the  Lagrange  coordinate  system 
is  considered  as  three  separate  ducts  joined  through  appropriate  boundary  conditions  at  the  end  of  each  duct. 
The  finite-difference  grid  in  each  of  the  ducts  is  moving,  which  requires  the  boundary  zones  to  be  modified  at 
each  calculation  cycle.  For  example,  with  the  shock  wave  moving  down  duct  1 into  ducts  2 and  3,  as  shown  in 
Figure  4,  the  boundary  zone  of  duct  1 is  moving  into  the  junction,  and  the  boundary  zones  of  ducts  2 and  3 
are  moving  away  from  the  junction;  this  requires  the  mass  to  be  removed  from  the  duct  1 boundary  zone  and 
added  to  the  boundary  zones  of  ducts  2 and  3.  This  boundary  zone  modification  (changing  the  mass  of  the 
zone)  is  done  in  a manner  that  satisfies  conservation  of  mass.  In  addition  to  conservation  of  mass,  the  conserva- 
tion of  energy  and  the  flow  losses  through  the  junction  are  also  considered.  Conservation  of  energy  is  provided 
by  assuming  a quasi-steadv  state  flow  exists  which  allows  use  of  a constant  total  enthalpy;  this  is  the  same 
method  previously  described  that  is  used  in  determining  flow  through  a duct  inlet.  The  flow  losses  through 
the  junction  are  determined  from  experimental  data  (Reference  8). 

The  computation  procedure  used  at  the  T-junction  is  an  iterative  one.  A value  for  the  static  pressure 
at  the  exit  to  duct  1,  pj1  N in  Figure  4.  is  determined  during  each  calculation  cycle,  which  allows  satisfaction 
of  the  following  conservation  of  mass  relations; 


ft"  » Am,  ♦ Am,  * Am, 


16"!  < m?  x 10'7 


(22) 


where  Am,,  Anij,  and  Am,  arc  the  changes  in  the  mass  of  the  boundary  zones  of  ducts  1,  2,  and  3 during  the 
nlh  time  cycle,  respectively,  and  m"  is  the  total  mass  of  the  boundary  zone  of  duct  1.  Through  an  iteration 
procedure  the  deviation  of  ft™  from  zero  is  determined  to  within  a small  bound.  The  following  auxiliary 
relationships,  in  finite-difference  form,  are  used  to  determine  Am,. 


Am,  * 

-ai  pi^n-i/z  (*?:4  ■ li  J 

(23a) 

n*l 

*I,N 

= *n  * Ar"*1  i.n+,/2 

XI.N  * UI.N 

(23b) 

n*l/2 

UI.N 

_ „n  . A,n»l/2_n 

* U1,N  ♦ *I.N 

(23c) 

»I.N  - 

2A , A"-/a  ‘ P?  N) 

(23d) 
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flow  direction 


Figure  4.  T-junction  finite-difference  grid  notation. 
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^l,N-l/2 


vl,N-l/2 


(23f) 


In  the  above,  L,  is  the  length  of  due:  1,  A,  is  the  cross-sectional  area  of  the  bi.:in  jary  zone,  and  other 
symbols  are  as  previously  defined.  The  acceleration  cf  the  Nfh  interface  of  duct  1 at  the  n1*1  calculation 
cycle  (or  nlh  time  step)  is  represented  by  aJ  N. 


The  inlet  pm  turn  to  duct*  2 and  3,  p2  , and  p " , , mpectively,  are  determined  from  the  exit 
pretsure  of  duct  1 through  the  data  of  Reference  8.  Losses  are  accounted  for  in  the  same  manner  as 
previously  discussed  for  the  inlet  simulation  of  duct  I,  i.e.,  through  change  in  static  pressure  following 
a normal  shock  wave  propagating  through  the  junction.  The  pressures  pj,  j and  p".  t depend  upon  p"  ^ 
through  relationships  similar  to  Kquation  18  but  with  different  empirical  constants.  (In  the  computer 
program  these  relationships  are  input  in  table  form.)  The  mass  transferred  from  duct  1 to  duct  2.  Am2, 


is  determined  in  terms  of  p2i  t (and,  therefore,  in  terms  of  p"  N implicitly)  by  the  relationships  given  below. 
The  mass  transferred  to  duct  3,  Antj,  is  determined  by  relationships  of  identical  form  but  with  pj_ , in  place 
of  p2  | and  index  2 replaced  by  3 in  other  appropriate  symbols. 

Amj  • A*  Pi*hi  x2r»* 

(24a) 

*Vi  m *5.1  ♦ 

(24b) 

n*  1 H _ „it  . »,n*l/2.n 

“2.1  " “2.1  * at  *2.1 

(24c) 

« . . /P2.3/2  * P2.l\ 

#2  ‘ ‘ 2*\  m2  J 

(24d) 

/wn*l  „n*l\ 

v-i  . r12  **•*  u 

\ „5  / ’ 

(24e) 

■ rnr 

v2.J/2 

(240 

In  the  above.  A,  is  the  cross-sectional  area  of  the  boundary  zone  of  duct  2.  x2  2 is  determined  using 
Equation  7b,  and  Pj  j/j  represents  p2  ♦ q2  j/2  Utilizing  Equations  23,  24,  and  the  equivalent  of 

Equations  24a  through  24f  for  Am  j,  Equation  22  ir  satisfied  by  iteration  of  the  pressure  pjN  for  each 
calculation  cycle. 

Since  mass  is  added  to  the  boundary  zones  of  ducts  2 and  3 from  duct  1 , the  mass  of  the  three 
zones  involved  is  adjusted  at  each  calculation  cycle.  In  addition,  the  internal  energy  and  static  pressure  of 
the  boundary  zones  of  ducts  2 and  3 arc  adjusted  to  account  for  effects  of  the  internal  energy  of  the  added 
mass.  Assuming  quasi-steadv  state  flow,  and  therefore,  constant  total  enthalpy,  and  using  a mass  weighted 
averaging  method  yields  the  following  relationships  for  the  adjusted  variables. 

„n*  1 _ „n+ 1 

*l,N-l/2  M.N-l/2 

(25a) 

_n*  1 • 

un*l  . Pl,N-l/2  1 /__n*  t \ 

htl  * Si.N-I/2  ♦ , + 2 p.*; 

Pl.N-l/2 

(25b) 
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(25c) 


||+  t * 


(>2.3/2 
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m"* 1 * itij  ♦ Am, 
mj*1  " n»2  ♦ Am2 


mj  +■  Amj 


£2*3/2  ’ (m2  e2,*3/2  " Am2  ^*2)  “^77 

ss:s'/i  * (,niei*i/i  ~ Am3Aej) 

\ / mi 


JJt.N-l/2 


_n*  1 

Pl.N-l/2 


tfUn  * <7  - ^^nPz^n  (25k> 

JJSN/*  * <*  - "GAnKln  (25,) 

The  barred  quantities  are  adjusted  variables,  and  values  for  P,*4-l/2’  c?.N-l/2’  P2.3/21  c2.3/2’  P3.3/2’  *n^ 
c»-j'/2  are  determined  from  Equations  7e  and  7f.  Values  for  other  variables  arc  obtained  from  Equations  23  and  >4. 

The  mass  of  the  boundary  zone  of  duet  1 increases  and  the  mass  of  the  boundary  zones  of  ducts  2 
and  3 decrease.  Therefore,  a limit  has  to  be  set  on  the  change  in  the  mass  of  these  zones.  A rezoning  method 
is  used  that  merges  the  boundary  zone  of  duct  1 with  the  adjacent  zone  when  the  boundary  zone  mass 
becomes  one  half  of  its  original  (cycle  I)  value.  The  boundary  zones  of  ducts  2 and  3 arc  divided  into  two 
zones  of  equal  mass  when  their  mass  becomes  double  the  original  (cycle  1 ) value.  This  rezoning  method  will 
be  discussed  in  more  detail  in  the  section  on  duct  rezoning. 

Equations  23  through  25  are  restricted  in  use  to  the  flow  and  shock  wave  directions  shown  in 
f igure  4.  i.e..  flow  from  duct  1 into  ducts  2 and  3.  This  is  the  initial  flow  pattern  for  a shock  propagating 
into  the  system  inlet.  Shock  reflections  from  the  debris  pit  end  and/or  the  blast  valve  (Figure  2)  can  cause 
the  flow  directions  to  change.  The  computer  code  subroutine  is  written  for  the  case  given  in  detail  above, 
and,  therefore,  computations  will  terminate  when  a reflected  shock  wave  arrives  at  the  T-junction.  Subrou- 
tines can,  however,  easily  be  a ided  to  the  computer  code  to  compute  propagation  of  reflected  shock  waves 
through  this  junction  since  tin  computation  method  is  similar  to  the  case  given  here;  only  the  directions  of 
mass  transfer  and  flow  loss  relationships  change.  The  main  limitation  this  restriction  places  on  the  use  of  the 
computer  code  in  its  present  form  is  the  ability  to  compute  long  er  ugh  in  time  to  allow  flow  reversal  in  both 
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the  duct  to  the  debris  pit  and  the  duct  to  the  blast  valve.  For  example,  if  the  blast  valve  duct  is  long  compared 
to  the  debris  pit  duct,  a reflected  shock  wave  from  the  debris  pit  could  arrive  at  the  T-junction  before  the  shock 
wave  in  the  blast  valve  duct  arrives  at  the  valve.  This  limitation  is  overcome  by  computing  two  geometries,  the 
desired  geometry  to  obtain  the  debris  pit  solution  and  a second  geometry  that  replaces  the  debris  duct  with  a 
very  long  constant  area  duct  to  obtain  the  blast  valve  duct  solution.  This  method  is  relatively  accurate  since 
the  debris  pit  flow  docs  not  significantly  affect  the  flow  in  the  blast  valve  duct  until  a reflected  shock  arrives 
at  the  T-junction  (see  Kesults  Section).  In  this  manner  the  primary  problem  is  solved  which  is  the  temperature 
and  pressure  distribution  in  the  ducting  system  up  to  the  time  flow  reversal  occurs  at  the  blast  valve. 

Duct  Exit  Boundary  Condition*.  The  boundary  condition  at  the  exit  end  of  ducts  2 and  3 is  considered 
a closed-end  condition.  In  finite-difference  form  the  boundary  conditions  are  written, 


■l.N 


<26«> 


a" 

*J.N 


0 (26b) 

where  a£  N and  a”N  are  the  accelerations  of  the  N,h  or  last  interface  of  ducts  2 and  3,  respectively. 

The  Finite-Difference  Grid  Keznning  Methods 

Changing  the  finite-difference  grid  in  each  duct  by  adding  or  combining  zones  and  renumbering  the 
grid  system-denoted  rezoning— is  required  periodically  during  a calculation  to  allow  mass  flow  into  or  out  of 
a duct.  For  the  configuration  of  Figure  2 rezoning  is  required  for  calculating  mass  flow  into  duct  I and  also 
through  the  T-junction.  Another  type  of  rezoning  that  is  employed  in  the  Ugrange  computer  code  is  the 
removal  of  excessively  compressed  zones  by  combining  these  zones  with  an  adjacent  zone.  This  is  necessary 
to  control  the  minimum  size  of  the  time  step,  which  is  determined  for  each  calculation  cycle  by  the  smallest 
zone. 

Inlet  Duct  Rezoning.  At  the  inlet  of  duct  1 a new  boundary  rone  is  established  whenever  the  first 
interface  has  moved  a distance  equivalent  to  the  mass  of  the  initial  boundary  zone.  To  establish  this  zone 
the  state  of  the  gas  at  the  zone  center,  the  velocity  of  the  zone  interface,  and  the  zone  size  have  to  be  specified. 
The  zone  mass  is  determined  from  the  gas  state  and  zone  size,  and  the  acceleration  of  the  interface  is  deter- 
mined from  the  zone  pressure  with  adjacent  zone  pressures.  The  specification  of  initial  values  for  these 
quantities  is  necessarily  approximate;  however,  test  calculations  with  the  computer  code  show  that  rezoning 
causes  the  shock  pressure  to  be  within  5%  of  the  correct  value.  The  computed  shock  pressure  is  always  above 
the  correct  value,  yielding  a conservative  answer. 

The  basic  assumptions  that  are  made  to  establish  the  initial  vrlucs  for  the  new  zone  parameters  are, 
referring  to  Figure  3,  aj  = 0 and  JJij/j  * Neglecting  the  friction  term  (which  is  small)  in  the  first  part 

of  Equation  7a  and  using  the  following  relationship  for  the  acceleration  of  the  boundary  interface, 


*1 


yields  the  relationship 


(27) 


JTj/2  * Pj/2 


(28a) 
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<2«b> 


Again,  Pj/2  * p"/2  ♦ q"/2-  The  altovc- mentioned  assumptions,  therefore,  imply  that  the  acceleration  of 
the  boundary  interface  docs  not  change  when  a new  boundary  /one  is  created.  Of  several  rezoning  methods 
tested  the  above  yielded  the  most  accurate  results. 

The  criterion  for  establishment  of  the  new  boundary  /one  is 


*o  > *» 


(29) 


where 


(30) 


<31 ) 


and  hrj  is  given  by  Kquation  14. 

The  initial  values  of  parameters  for  the  new  boundary  zone  of  duct  I based  on  Kquations  28  are  as 

follows: 
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JP3/2  * P3/2 

(32d) 

£°m  “ 
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„ <7  ” 1 >£3/2 
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A,  xj 
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When  initial  values  have  been  established  for  the  new  zone  parameters,  all  the  zones  of  the  finite-difference 
grid  are  renumbered  to  include  this  additional  zone. 

T-Junction  Rczoning.  Since  mass  is  flowing  through  the  T-junction,  the  boundary  zones  at  this 
junction  change  mass  with  each  calculation  cycle.  The  exit  boundary  zone  of  duct  t decreases  in  mas*,  and 
the  entrance  boundary  zones  of  ducts  2 and  3 increase  in  mass.  This  change  in  zone  mass  eventually  requires 


14 


(hat  the  boundary  zones  be  altered.  The  procedure  that  is  used  in  the  computer  code  is  similar  to  the  duct  I 
inlet  rezoning  described  above,  i.e.,  a new  zone  is  established  periodically,  and  the  finite-difference  grid  renum- 
bered accordingly.  The  exit  boundary  zone  of  duct  1 is  combined  with  its  adjacent  zone  to  form  a new 
boundaiy  zone  when  its  mass  becomes  one  half  of  its  original  mass.  The  inlet  boundary  znr.es  of  ducts  2 
and  3 are  split  to  form  two  zones  when  their  mass  becomes  double  the  original  mass. 

The  criterion  for  establishment  of  a new  exit  zone  for  duct  1 is 


"St-l/J 


< 


1 

7 


mN-l/2 


(33) 


where  mj5,_ ( /a  represents  the  zone  mass  at  time  zero.  The  relationships  for  establishing  values  of  parameters 
for  the  new  zone  are  given  below. 
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In  the  above.  ejJ_1/2,  etc.,  represents  the  new  zone  value. 

The  criterion  for  establishment  of  a new  entrance  zone  for  duct  2 or  3 is 

iHj/2  ^ (35) 

where  m j/2  represents  the  zero  time  value.  The  relationships  for  establishing  the  two  new  entrance  zones  are 
as  follows. 
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These  relationships  apply  to  establishing  new  zones  at  the  entrance  of  both  duets  2 and  3. 

Rezoning  to  Reduce  Running  Time.  To  aid  in  reducing  computer  running  time  it  is  expedient  to 
avoid  zones  becoming  excessively  compressed.  The  calculation  time  step  is  based  upon  the  smallest  or  most 
highly  compressed  zone.  When  the  time  step  becomes  too  small,  it  is.  therefore,  necessary  to  remove  the 
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iiMJIett  zon*  (the  zone  upon  which  the  time  step  is  bucd)  by  combining  this  zone  with  an  adjacent  zone.  A 
criterion  for  removal  of  the  small  zone  is  established,  and  relationships  for  combining  this  zone  with  an  adjacent 
zone  are  formulated. 

The  criterion  for  removal  of  an  excessively  compressed  zone  is 

At"* I71  < j At®  (37) 

where  At®  is  the  initial  or  zero  time  value  of  At"*  ,/2.  The  relationships  to  establish  a new  zone  comprised 
of  the  zone  being  removed  and  an  adjacent  zone,  with  x"  the  location  of  interface  between  the  zones,  are: 
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The  quantities jj",  uj\  etc.,  represent  values  of  parameters  for  the  newly  created  zone,  and  A"  is  the  cross-sectional 
area  at  the  interface  of  the  zones  being  combined.  A zone  is  removed  whenever  the  minimum  time  step,  At"* 1/2, 
based  on  that  zone  becomes  smal'er  than  one  half  the  value  of  At"* 1/2  at  time  zero  as  expressed  by  Equation  37. 

The  Equation  of  State  Options 

Two  equations  of  state  subroutines  are  available  in  the  Lagrange  computer  code:  (a)  the  equation  of 
state  for  an  ideal  gas  that  is  accurate  for  temperatures  up  to  1,000°K,  and  (b)  a real  gas  equation  of  state  that 
is  accurate  for  temperatures  up  to  24,000°K.  The  basic  theory  upon  which  these  subroutines  are  based  is  given 
below. 
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Th*  Perfect  Cm  Equation  of  State.  The  equation  of  state  for  an  ideal  gas  ean  be  written  in  the  form: 


P 


(39) 


where  p is  the  pressure,  K ihe  particular  gas  constant,  T the  absolute  temperature,  and  V the  specific  volume. 
Itccjuse  of  the  computational  advantages  offered  by  the  use  of  Equation  3V,  most  problems  of  shock  trans- 
mission through  gases  are  solved  using  the  ideal  gas  law.  Moreover,  the  ideal  gas  law  yields  simple  analytical 
relationships  between  the  various  shock  parameters.  The  thermally  and  caloric,  llv  perfect  gas  is,  however, 
an  idealization,  and  real  gases,  depending  on  their  temperature  and  pressure,  will  deviate  from  it  to  a 
varying  degree.  It  has  been  noted  from  experience  < References  12  and  13)  that  the  ideal  gas  law  predicts 
the  shock  flow  parameters  with  good  accuracy  up  to  temperatures  of  However,  at  higher  temper- 

atures, the  number  ol  particles  per  unit  mass  and,  hence,  the  average  molecular  weight  of  the  gas  may  change 
due  to  molecular  dissociation,  chemical  reactions,  and  ionization.  Thus,  when  computing  flow  parameters 
tor  shocks  through  gases  at  high  temperatures  and  moderately  high  pressures,  the  ideal  gas  equation  must  be 
modified  to  include  the  contributions  by  additional  degrees  of  freedom  to  the  energy  of  the  gas  molecules 
that  were  not  excited  at  low  temperatures. 

The  Real  <ias  Equation  of  State.  The  translational  and  rotational  degrees  of  freedom  of  molecules 
are  excited  at  very  low  temperatures  up  to  the  order  of  IO°K,  whereas  the  molecular  vibrations  absent  in  mono* 
tonic  gases  are  fulls  excited  at  much  higher  temperatures  of  'he  order  of  2.2  for  Oj  and  3,340°K  for  Nj. 

I he  contribution  by  the  vibrational  excitation  at  lower  temperatures  can  be  computed  by  the  method  of 
partitions  described  in  Reference  1 2.  t hus,  for  a diatonic  gas  with  no  molecular  dissociation  and  ionization 
present,  ■>.  the  adiabatic  exponent,  is  7/5  with  unexcited  vibrational  mode,  while  with  fully  excited  molecular 
vibrations  > = 9/7.  l-or  a monotonic  gas  > holds  a constant  value  of  5/3. 

It  should  be  staled  here  that  the  temperature  range  for  which  the  molecular  vibrations  stay  fully, 
excited  is  not  very  w ide.  This  is  In-cause  the  molecular  dissociation  and  chemical  reactions  frequently  begin 
at  temperatures  at  which  the  molecular  vibrations  contribution  to  internal  energy  reaches  its  classical  limiting 
value.  At  temperatures  of  the  order  of  several  thousand  degrees  < 3.(MWt°K).  the  molecules  of  diatonic  gases 
dissociate  into  atoms.  The  process  of  dissociation  of  a molecule  requires  a large  amount  of  energy,  thus 
affecting  the  thermods  mimic  properties  of  the  gases  appreciably.  At  standard  atmospheric  density,  molecular 
dissociation  in  air  is  noticeable  ,u  2.9KI»‘,K  due  primarily  to  breaking  up  of  the  IN  molecule.  The  internal 
energy  of  a partially  dissociated  gas  is  the  sum  of  the  energy  of  the  umlissociatcd  molecules,  the  energy 
contained  in  dissociated  atoms,  and  the  energy  required  to  dissociate  unexcited  molecules.  Because  of  the 
high  energy  required  to  dissociate  a molecule  into  its  atoms,  the  energy  of  dissociated  gas,  even  for  small 
degrees  ol  dissociation,  is  very  much  greater  th  in  that  in  the  absence  of  dissociation,  l-or  monotonic  gases 
effects  due  to  molecular  dissociation  arc  absent. 

Another  factor  that  changes  the  thermodynamic  ptopcrtics  of  a mixture  of  gases,  such  as  air  at  high 
temperatures,  is  the  occurrence  of  chemical  reactions  between  its  constituents.  Air  at  temperatures  above 
l.5tMl"K  undergoes  oxidation  of  a portion  of  its  nitrogen  to  form  nitric  oxide.  This  is  a reversible  reaction 
and  requires  a high  activation  energy  (21.4  kcal/molc).  The  time  to  attain  equilibrium  is  very  high  at  temper- 
atures below  1 .5<>o"K:  however,  at  temperatures  of  the  order  of  3,<HH>”K  and  above  the  chemical  equilibrium 
is  reached  in  10"*  seconds  or  less.  In  other  words,  the  chemical  composition  of  air  around  3.(H)0°K  and  above 
contains  nitric  oxide  as  one  of  its  components,  which  changes  the  thermodynamic  properties  of  air  somewhat. 

Next,  ionization  of  atoms  or  molecules  of  gases  in  air  begins  around  K.OOO’K.  The  degree  of  ioniza- 
tion increases  with  temperature.  When  the  air  temperature  is  of  the  order  of  3(),000°K.  all  of  its  atoms  are 
singly  ionized.  At  temperatures  above  3<*.0<><>”K.  the  air  undergoes  second  followed  by  third  ionization 
processes.  Since  this  investigation  deals  with  air  temperatures  of  24,<HIO('K  and  below,  only  first  ionization 
of  air  particles  will  he  considered  here.  Unlike  molecular  dissociation,  the  internal  energy  of  the  ionized  gas 
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is  a sunt  of  the  thermal  energy  of  its  particles  (atoms,  ions,  and  electrons)  and  the  energy  tequired  to  remove 
electrons  from  the  atoms  and  ions.  The  internal  energy  of  a partially  ionized  gas,  therefore,  will  be  the  energy 
of  the  ionized  and  unionized  gas  particles. 

To  examine  the  behavior  of  thermodynamic  properties  of  air  with  temperature,  three  of  its  primary 
properties,  i.e.,  the  dimensionless  internal  energy  (e/R  T),  enthralpy  (h/R  T),  and  the  adiabatic  exponent  (y), 
were  obtained  from  the  data  given  in  Reference  13.  These  parameters  were  obtained  between  a temperature 
range  of  300  and  24,000°K.  A tabulation  of  the  data  at  two  different  densities,  namely,  at  the  standard  air 
density  (p0  * O.OOI293  gm/cm3)  and  at  1/10  of  the  standard  density,  is  given  in  Table  1.  The  variation  of 
the  parameters  e/R  T and  h/R  T with  temperature  is  shown  by  the  curves  of  Figure  5.  The  temperature  lines 
at  which  the  contributions  by  various  degrees  of  freedom  become  accountable  arc  dearly  marked  on  the 
curves.  It  is  obvious  from  Table  1 and  from  the  curves  that  the  energy  functions  e/R  T and  h/R  T undergo 
a wide  variation  as  the  air  temperature  increases  from  300  to  24,<tOO°K.  Correspondingly,  the  adiabatic 
exponent,  y,  goes  through  a wide  variation  from  the  classical  value  of  1.4  at  300°K  to  around  1.15  at 
lO.OwP'K.  Figure  6 shows  plots  of  y versus  temperature  at  the  t»v«*  selected  air  densities. 


Table  1.  Thermodynamic  Properties  of  Air  (After  Reference  1 3) 


Thermodynamic  Properties"1 

Temperature,  T 
(°K) 

At  Standard 
Density,  po* 

At  1/10  of  Standard 
Density.  P„ 

c/RT 

h/RT 

7 

e/RT 

h/RT 

y 

300 

2.5 

3.5 

1.4 

2.5 

3.5 

1.4000 

1,000 

2.630 

3.63 

1.38 

2.63 

3.63 

1.3800 

2,000 

2.962 

3.963 

1.3378 

2.966 

3.967 

1.3375 

3,000 

3.395 

4.402 

1.2967 

3.677 

4.699 

1.2780 

4,000 

4.345 

5.408 

1.2447 

5.311 

6.445 

1.2135 

5.000 

5.138 

6.283 

1.2290 

5.746 

6.945 

1.2087 

6,000 

5.502 

6.700 

1.2178 

6.517 

7.784 

1.1944 

7, otto 

6.623 

7.533 

1.2028 

8.732 

10.163 

1.1639 

8,000 

7.744 

9.144 

1.1808 

11.67 

13.349 

1.1439 

1 2,000 

12.18 

14.139 

1.1608 

13.47 

15.531 

1.1530 

18, (KM) 

13.04 

15.323 

1.1751 

18.13 

20.906 

1.1531 

24,000 

16.09 

18.96 

1.1781 

21.83 

25.432 

1.1650 

a e is  the  specific  internal  energy,  and  h is  the  specific  enthalpy  of  air. 

The  adiabatic  exponent  y ■ h/e. 

* p„  * 1.293  x 10'3  gm/cm3. 
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Figure  6.  Variation  of  the  adiabatic  exponent  with  temperature  for  air. 

where  Z(p,T)  » f(a,a')  (4!) 

a is  the  degree  of  molecular  dissociation,  and  a 1 is  the  degree  of  particle  ionization.  The  quantity  Z,  commonly 
known  as  the  “compressibility  factor,"  is  a function  of  pressure  and  temperature  of  the  gas.  At  temperatures  of 
1,(MH)°K  or  less  and  moderate  pressures  the  perfect  gas  equation  gives  satisfactory  results.  Thus,  Z * 1.0  at 
these  temperatures  and  pressures,  yielding  the  perfect  gas  equation.  The  quantity  Z for  air  at  a given  temper- 
ature and  density  is  found  from  standard  tables  given  in  Referenees  1 3 and  1 5.  These  are  included  as  Table  2 
for  easy  reference. 

The  Rankine-Hugoniot  Relationships  for  High  Temperature  Shock  Waves.  When  moderately  high 
overpressures  and  high  temperature  shocks  propagate  through  cold  air,  the  thermodynamic  properties,  such 
as  the  internal  energy  function  (c/R  T),  the  enthalpy  function  (h/RT),  and  the  adiabatic  exponent  <7  = h/c), 
vary  nonuniformly  with  temperature  (Figures  5 and  6).  The  nonuniform  variation,  as  discussed  earlier,  is 
attributed  mainly  to  molecular  dissociation  and  particle  ionization.  Thus,  the  Rankine-Hugoniot  relationships, 
because  of  their  dependence  upon  the  above  parameters,  change  in  a nonuniform  manner.  Therefore,  it  is  not 
possible  to  define  simple  analytical  expressions  describing  the  pressure,  density,  and  temperature  ratios  across 
the  shock.  There  are  some  iterative  procedures,  which  are  based  upon  the  tabulated  value  of  the  compressibility 
factor  Z and  the  adiabatic  exponent  7,  which  render  semi-anaivtical  expressions  for  the  shock  wave  parameters 
112.13,141. 

Gilmore  1 1 31  plots  the  density  and  temperature  ratios  (p2/p,  and  T2/T, ) across  the  shock  in  the 
standard  llugoniot  form,  i.e.,  against  the  pressure  ratio  (P2/P, ) for  density  values  behind  the  shock  p2  of  10, 

1.  10  *,  10'2,  10  J,  and  10  of  the  standard  air  density  po  = 0.001293  gm/em2.  Gilmore  in  his  plots  defined 
suffix  1 to  refer  to  the  quantities  ahead  of  the  shock,  and  suffix  2 was  used  for  quantities  behind  the  shock. 


The  plot*  are  based  upon  a normal  shock  wave  propagating  into  cold  air  at  a temperature  Tj  of  273. 2°K.  The 
plots  alto  show  the  variation  of  the  Mach  number  of  shock  as 


u^  /(p2/pl)(p2/p,  - I) 

cl  V Ti<P2/P|  - I) 


(42) 


with  the  pressure  ratio  p2/p|.  The  quantity  u,  is  the  speed  of  the  shock,  and  ct  is  the  speed  of  sound  in  the 
cold  air.  These  Kankinedlugoniot  plots,  which  include  the  real  gas  effects,  arc  given  in  Figure  7 for  normal 
shocks  traveling  through  air  at  a density  of  0.001293  gm/cm J. 

Computations  Using  the  Real  Gas  Equation  of  State.  The  l.agrange  computer  program  has  ».t 
equation  of  state  algorithm  built  into  it.  The  algorithm  based  upon  a value  of  the  internal  energy,  e,  and 
the  air  density,  p,  computes  the  adiabatic  exponent,  y.  from  the  tabulated  values  of  Gilmore  |I3| . From 
known  values  of  e,  p,  and  y the  value  for  the  pressure  is  computed  using  the  relationship 


p ■ pe<7  - 1)  (43) 

From  the  known  value  of  p,  the  algorithm  refines  the  previously  calculated  values  of  c and  p for  further 
computing. 

The  algorithm  has  a built  in  numerical  routine  to  obtain  a correct  value  of  air  temperature.  The 
method  is  based  upon  computing  the  air  temperature  from  known  p and  p using  Equation  40.  Since  the 
quantity  Z is  a function  of  pressure  p and  temperature  T,  calculation  of  temperature  using  Equation  40  is, 
therefore,  an  iterative  process.  The  numerical  routine  performs  the  required  iteration  accurately,  using 
the  tabulated  values  (Table  2)  of  Z given  by  Gilmore  1 1 3 1 . 


Figure  7.  Rankine-llugoniot  curves  including  real  gas  effects  for  shotk  waves  in  air. 
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Table  2.  Value*  of  the  ('(impressibility  Factor  Z(p,T)  (After  Reference  U) 


Temperature,  T 
(°K) 

(.'(impressibility  Factor  for  Density  Ratio,  p/p„. 

•f- 

to 

10° 

10* 

HI'2 

Hf* 

1(T4 

10  5 

1()'6 

1,000 

1.0000 

1.(8881 

1.0000 

1.(8881 

1.(888) 

1.0000 

1.(8881 

1.(888) 

HZ « 

1.0000 

1.(888) 

1.(8811 

1.(8817 

1.(8117 

I.OOJ5 

1.0170 

1.0487 

■ 

1.(8123 

1.(8172 

1.0610 

1.1341 

1.1912 

1.2069 

1.2109 

4,000 

1.0226 

1.0628 

1.1896 

1.2282 

1.2715 

1.3883 

5, Otto 

1.0728 

1.1448 

1.1990 

■ 

1.4835 

1.7752 

1.9544 

6,000 

1.1275 

1.1984 

1.2669 

1.4(811 

■ 

1.9143 

1.9864 

2.0051 

7,(881 

1.1761 

1.2699 

1.4315 

1.7162 

1.9970 

2.0288 

2.1105 

8,000 

1.2390 

1.4001 

1.6794 

1.9259 

2.(8110 

2.0487 

2.1702 

2.5(817 

12.000 

1.708 

1.957 

2.062 

2.224 

2.636 

3.365 

3.865 

KSSB 

18,000 

2.063 

2.284 

2.775 

3.518 

3.910 

3.984 

4.078 

KSH 

24,(NNt 

2.322 

2.865 

3.602 

3.930 

4.143 

4.778 

5.623 

I n 

Tcit  Calculation*.  To  check  the  accuracy  of  the  algorithm,  four  trial  runs,  using  the  Lagrange  computer 
program,  for  shocks  traveling  into  still  air  at  standard  density  and  temperature  through  a straight  duct  were 
carried  out  at  pressure  ratios  of  50:55.  101:77,  502:45.  and  2008:8  across  the  shock.  The  computed  shock 
Mach  number.  Mv  and  the  density  ansi  temperature  ratios,  p2lp , and  T2/T,,  across  the  shock  are  given  in 
Table  3.  These  values  arc  also  plotted  on  the  Kankinc-llugoniot  curves  of  Figure  7 for  comparison  with 
Gilmore's  accurate  values.  It  can  be  seen  from  Figure  7 that  the  shock  parameters,  namely  the  Mach  number, 
and  the  density  and  temperature  ratios  across  it  arc  computed  with  good  accuracy  by  the  Lagrange  computer 
program.  For  completeness,  the  values  of  the  adiabatic  exponent  y and  the  compressibility  factor  Z arc  also 
listed  in  Table  3. 


Tabic  3.  Shock  Parameters  Computed  Using  the  Real  Gas  Liquation  of  State11 


Pressure 
Behind 
Shock,  pj 
(psia) 

Pressure 

Ratio. 

lypi 

Shock 
Mach 
Number, 
M,  • u,/c, 

Density 
Behind 
Shock, 
(gm/um* ) 

Density 

Ratio. 

p2ipi 

Temperature 

Behind 

Shock. 

t2  ("Kl 

Temperature 

Ratio. 

t2/t, 

Adiabatic 
ttxpo.ient,  y 

Compressibility 
Factor,  7. 

743.1 

1 

mm 

7.447  x III  * 

6.079 

mm 

mm 

1.33 

1.00 

1.496.0 

mum 

9.051  x HI'3 

7.389 

mSSSm 

1.28 

1.02 

7.3K6.0 

1.341  x 111'2 

1095 

9.351 

ilSSSi 

1.19 

1.41 

29.530.0 

2.00M.H 

; 

1.399x  Uf2 

11.42 

22.470 

■fil 

1.19 

2.26 

Pressure  (pj).  temperature  (Tj),  and  density  <p  j)  for  the  undisturbed  air  are  taken  as  14.7  psia. 
2KM.2°K,  and  1.225  x 10*^  gm/cm*.  respectively. 
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RULKKIAN  COMPUTER  CODE  DESCRIPTION 


The  Buie  Equation* 


The  Kulcrian  code  computes  the  flow  variables  in  a duct  through  which  a normal  shock  wave  is 
passing.  The  flow  variables  are  computed  as  functions  of  time  and  position.  The  shock  wave  is  considered 
to  be  transmitted  into  the  duct  by  a blast  wave  passing  over  the  duct  entranci  The  computer  code  allows 
the  shock  wave  to  traverse  the  duct,  be  reflected  from  the  closed  end  of  the  duct,  and  then  re-traverse  the 
duct  in  the  opposite  direction. 

The  conservation  equations  are  written  in  the  vector  form. 


3U 

dt 


3F(U) 
* flx 


♦ G(U) 


0 


(44) 


where  vectors  U,  K(U).  and  G(U)  are  given  by 


The  quantities  M,  N,  and  E are  defined  as  follows: 


(45«) 


(45t>) 


(45c) 


M * p u A 
N * pA 
E * pie  + 1/2  u2) 
t = time 
x » axial  distance 
p - density 

c = internal  energy  per  unit  mass 

p » pressure 

= p(e,p) 

= p eiy  - I ),  for  a polytropic  gas 
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f • wall  friction  coefficient 


■ 4rw/(l/2pux) 

u » velocity 

rw  * shear  stress  .it  duct  wall 
A « duct  cross-sectional  area 
II  * duct  hydraulic  radius 

Substitution  ot  the  vectors  U,  F(U),  and  G(U)  into  Equation  44  yields  the  equations  for  the  conservation  of 
mass,  momentum,  ami  enemy.  The  gas  in  the  duct  is  considered  to  be  inviscid  (except  at  the  duct  wall)  and 
thermally  nonconducting.  Pressure,  density,  and  internal  energy  are  related  through  the  perfect  gas  equation 
of  state. 

Because  the  flow  in  the  duct  is  time  dependent,  steady-state  values  of  the  wall  friction  coefficient,  f, 
cannot  properly  Ik-  used.  In  the  Kulerian  code,  the  wall  friction  coefficient  is  computed  from 

f • SKRc)s2  (4A) 

where  K,  is  the  Reynolds  number,  based  <m  duct  equivalent  diameter,  and  SI  and  S2  arc  empirical  constants. 
The  Finite-Difference  Equations 

A two-step  integration  technique,  that  of  I. ax  and  Wendroff  (Reference  2),  is  used  to  solve  the 
conscrvat-ot  equations.  Equations  44  and  45.  Referring  to  the  sketch  below,  all  the  properties  at  time  step 
n arc  known.  The  properties  at  time  step  n+1  and  location  j arc  to  be  computed. 


Axial  Distance  Along  Duct 
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First,  the  properties  at  two  intermediate  steps  {J— 1/2,  n+1/2)  and  (j+1/2,  n+l/2)  are  computed  from: 


u?-V/i'  *u")  - 4 (s  * <:)  (-H  - l) 

■ 4 • ur)  * t (cr..  * '■/) 

ur-V/,'  ■ T (op.,  • u,")  -■!(£-  «-,)(►■;  - ►:,) 

- 4 «;  fr  - uh)  * f (<:f-.  * <?) 


where  Fi,  Gj,  etc.,  denote  F(Uj),  G(Uj),  etc.,  and  gp  gj,  gj,  gj  are  pseudo-viscosity  factors,  defined  in 
Fquation  49.  With  the  properties  at  the  two  intermediate  steps  known,  the  properties  at  the  new  step 
(j,  n+l)  are  computed  from 

u?"  • u”  - 3"  - TV/?) 

* 4 (=)*;(“)•■  - 

, h-n*\n  . rn*l/2'i 

♦ Y \°j+l/2  + °n-l/2/ 

where  gj|  and  g|i  are  pseudo-viscosity  factors,  defined  in  Kquation  49.  The  pseudo-viscosity  factors  were 
obtained  following  the  generalized  procedure  of  Reference  4.  This  procedure  yielded  the  following 
equations: 

gj  » a,  [(* t/2y  - (cj*.  1/2)  2 J * a2  Uj"  I >2  (U)>  1 n ~ cj"  I ll) 

* ai  %"♦  1 12  1 12  * cj"  I/2J 

g|  * ■ j^°l  1/2  + ®2  ^2  U",  l/2  - Cj"  l/2J  ♦ Oj  ^2  U"*1/2  + Cj"  ,/2jJ  ( 


* 2 = “|  ♦ <*2  + a 3 


I ..  Iui"  t - uj"' 


o,  - --K, 


o2 ..  I Kz  ifo.-q)  + fo'-ci)' 


t50c) 


aS 


('"♦I/?) 


1 


ami  c > acuuitik  velocity. 

Values  of  Kt,  K2,  and  Kj  of  2.0  have  been  empirically  found  to  be  of  the  proper  magnitude.  The 
expressions  for  gg,  g'f , and  gj  are  obtained  by  replacing  j ty  j-l  in  Equations  49  and  50. 

The  stability  requirement  for  selecting  the  time  step  is: 


At 


Ax 

F(u  ♦ c) 


(51) 


where  F is  a factor  greater  than  unity.  Numerical  results  have  shown  that  too  small  a time  step  results  in 
oscillatory  solutions.  Thus,  the  value  of  F should  be  the  smallest  value  <>1 ) that  wilt  yield  a stable  solution. 
The  optimum  value  for  F was  found  to  be  1.5.  In  the  computer  code,  the  time  step  At  is  recomputed  for  each 
cycle,  using  the  minimum  value  of  lui  + c of  the  preceding  cycle. 

The  Boundary  Conditions 

The  pressure  at  the  entrance  of  the  duct  is  specified  as  a function  of  time  in  the  same  manner  as  in 
the  Lagrange  Code  (Equation  10): 

p<0,t)  - p.  ♦ [p<0.0>  - p.J  e (52) 


where  p<0,t)  » pressure  at  x * 0,  t * t The  initial  pressure  at  the  duct  entrance,  p(0,0),  is  computed  (as 
ir  the  Lagrange  Code,  Equation  18)  from: 

p(0,0>  - p,  ♦ 0.969(p,  - p,)0-804  (53) 

With  p<0,t)  known,  the  density  p(0,t)  is  computed  from  the  iscntropic  relationship 

p(0,,>  - P<0.0)[£M],/7  (54) 

where  7 is  the  ratio  of  specific  heats,  considered  constant  in  this  study. 

The  initial  value  of  density  p(0,0)  i»,  computed  from  the  Rankinc-Hugoniot  relationship 


P(0,0) 


Pj 


(7  ♦ 1) 


(7-1) 


p(0.0) 

P« 

p(0.0) 

P. 


(7-  1) 


(7+1) 


(55) 


where  p,  is  the  ambient  value  of  density.  Internal  energy  e is  computed  from  the  equation  of  state  (Equation 
43),  written  in  the  form 


e(0,t) 


p(0,t) 

p(0,t)(7-  1) 
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The  initial  velocity  u<0,0)  is  computed  from  the  moving  normal  shock  relationship 


u<0,0)  - 


where  us,  the  velocity  of  the  normal  shock  wave,  is  computed  from 


u 


rp<o.o) 

f-y  * t\  *■  Kwh 

p. 

IT  f ‘ * \ J * r 

co 

27 

(57) 


(58) 


where  c„  is  the  ambient  acoustic  velocity  Because  a rational  basis  for  specifying  u(0,  t)  a priori  is  not  known 
to  -lie  authors,  the  expediency  of  setting  u(rt,t)  equal  to  u(Ax,  t)  was  employed,  where  Ax  is  the  numerical 
incrcmcn'  of  x used  in  the  numerical  solution  of  the  conservation  equations. 

At  the  end  of  the  duct,  x = L.  the  rigid  wall  boundary  condition  can  be  imposed  as  u(L,  t)  * 0. 

Within  the  framework  of  the  numerical  integration  scheme,  this  was  accomplished  by  means  of  a virtual  station 
lievond  the  end  of  the  duct,  located  at  I.  + Ax.  By  setting  u(L  + Ax, t)  equal  to  -u<L- Ax, t),  the  proper 
boundary  condition  at  x * L was  obtained. 


NUMKKICAL  RESULTS 

l_i grange  Computer  Code 

The  results  from  the  l.agrangc  computer  code  are  presented  in  Reference  1 for  the  particular  case  of  a 
constant  area  duct.  A summary  of  these  results  is  reproduced  in  f igures  8 anil  9,  in  which  computed  data  are 
compared  with  the  shock  tube  experimental  data  of  Reference  <>.  These  comparisons  show  the  degree  of 
accuracy  of  the  computer  code  when  a constant  friction  factor  is  used  and  demonstrate  the  effect  of  wall 
friction  on  shock  wave  attenuation. 

Computed  results,  w hich  include  the  variable  area  option,  are  presented  in  l igure  10.  Shock  wave- 
pressure  change  through  a duct  cross-sectional  area  increase  is  given  for  outlet-to-mlet-area  ratios  up  to  10. 

I he  analytical  curve  is  compared  w ith  Stanford  Research  Institute  shock  tube  experimental  data  1 16 1 and 
with  data  from  the  joo-ton  I N I blast  experiments  of  KVTNT  1)1  Al.  PACK  1 5 1 . The  duet  diameter  at  the 
inlet  to  the  area  change  was  2 inches  in  the  SRI  experiments  and  24  inches  in  the  DIAL  PACK  experiments. 

The  comparison  demonstrates  the  adequacy  of  the  computer  code  to  predict  shock  strength  change  with  duct 
area  increase. 

Pressure-versus-time  waveforms  are  presented  in  figure  1 1 for  the  l.agrangc  Code  solution  with  a 
comparison  to  results  from  the  fulerian  Code  (discussed  in  a following  section).  The  reflected  wave  at  the 
closed  end  of  the  duct  is  also  shown  in  the  figure.  The  overpressure  of  this  reflected  wave  is  approximately 
It)".,  below  the  value  estimated  by  Kankinc-llugontot  relationships.  Studies  have  shown,  however,  that  the 
incident  wave  and  reflected  wave  overpressures  at  a rigid  Imundary  computed  for  a shock  wave  of  constant 
strength  by  the  l.agrangc  Code  agrees  well  with  Kankinc-llugontot  values:  this  implies  a small  reflected  wave 
inaccuracy  for  waves  with  rarefaction  behind  the  shock. 

The  computed  results  using  the  l.agrangc  Code  for  the  case  of  a l-Mt  nuclear  surface  wave  are  given  in 
f igure  I 2.  The  air  entrainment  system  is  a typical  configuration  except  for  possible  variations  in  the  scale.  The 
location  of  the  inlet  was  selected  as  1.500  feet  from  the  center  of  the  nuclear  hurst.  The  required  surface  wave 
parameters  were  obtained  from  Reference  7.  The  pressure  ami  temperature  values  show  n in  the  figure  depict 
the  state  of  the  air  in  the  system  at  a time  of  0.0108  second  after  arrival  of  the  nuclear  wave  at  the  inlet.  The 
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SRI  data,  Pj  * 33(1  paii.  Reference  16 
SRI  data,  pj  * 700  paia,  Reference  16 
DIAL  PACK  data,  pj  * 230  paia,  Reference  9 
Analysis,  moving  grid  solution,  f • 0,016 


■ I I I - I mi 

7 4 6 H 10 

Area  Ratio  (A./A, ) 

Figure  10.  Kffect  of  area  change  on  shock  pressure. 


t j » 1.10  msec 
D » 1,0  in. 
f • 0.016 


fixej  grid  solution  (Hulerian) 

1 

I moving  grid  solution  (l.agrangel 


l \ 


19.7  ft 


time  of  shock  arrival  at 
dosed  end,  x * 32.8  fi 


reflected  wave 

/Vj 


>1 


a 


26.2  ft  31.2  ft 


Figure  12.  Typical  computed  result*- Lagrange  code,  0.0108  second  after  arrival  of  shock 
wave  from  1-Mt  nuclear  explosion. 


shock  wave  and  cold-hot  gas  interface  (contact  surface)  locations  art  shown  as  a line  boundary.  The  computer 

program,  however,  spreads  the  shock  wave  over  approximately  five  finite-difference  grid  mesh  lengths,  due  to  i 

the  pseudo-viscosity  technique  used,  and  spreads  the  contact  surface  also  over  approximately  five  mesh  lengths. 

The  shocks  are  accurately  located  through  utilization  of  the  maximum  value  of  the  pseudo-viscosity  parameter,  j 

and  this  is  the  location  shown.  The  contact  surface  location  shown  is  the  location  of  the  zone  interface  across  . 

which  the  leading  downstream  temperature  jump  occurs  and  is  the  location  were  the  total  temperature  jump  ! 

will  occur  (neglecting  diffusion  effects)  in  the  actual  case.  Figure  12  is  an  example  of  the  computet  predictions  i j 

at  a particular  time  and  not  the  complete  results  of  the  computer  output.  For  example,  for  this  calculation  the  j 

time  history  of  all  parameters  was  available  at  several  locations  in  the  system.  These  time  histories  are  not 
shown  because  they  would  not  add  significantly  in  presenting  a sample  calculation.  The  computer  code  com- 
plete output  capabilities  are  explained  in  Appendix  A. 

When  a reflected  shock  wave  reaches  the  T-junction  (Figure  12),  the  computer  program  terminates  i 

computation  because  the  T-junction  programming  will  not  allow  reversal  of  the  direction  of  the  inlet  or  exit 
velocities  at  this  junction.  This  fact,  however,  does  not  seriously  limit  the  capabilities  of  the  computer  code. 

As  examination  of  Table  4 shows,  the  area  expansion  ratio  for  the  debris  pit  inlet  has  negligible  effect  on  the  1 

flow  parameters  in  the  inlet  duct  and  in  the  duct  to  the  blast  valve.  Therefore,  any  given  problem  can  be 

caicu'ated  in  two  steps  as  follows:  step  1,  calculate  the  inlet  duct  and  debris  pit  flows  rsing  the  configuration 

shown  in  Figure  lb;  step  2,  calculate  the  blast  valve  duct  flow  by  using  the  configuration  shown  in  Figure  2 j 

with  the  debris  pit  duct  and  debris  pit  replaced  by  a constant  area  duct.  This  constant  area  duct  should  be 
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sufficiently  long  to  allow  the  reflected  wave  in  the  blast  v.lve  duet  to  reach  the  T-junction  prior  to  the  time 
when  a reflected  wave  in  the  constant  area  duct  would  reach  the  T-junction.  This  allows  sufficient  completion 
of  the  primary  shock  propagation  problem  in  the  blast  valve  duct.  The  debris  pit  solution  (step  I),  however, 
will  lie  conservative  since  the  pressure  loss  through  the  T-junction,  which  is  approximately  20*K,  will  not  be 
included.  Neglected  in  the  ab  -"c  procedure  is  the  interaction  of  primary  and  reflected  shock  waves,  e.g..  the 
effects  of  j reflected  wave  from  the  debris  pit  propagating  into  and  down  the  blast  valve  duct.  Computation 
of  shock  wave  interaction  is  not  considered  important  since  flow  reversal  will  have  taken  place  in  all  ducts 
prior  to  any  shock  interactions  for  the  configuration  of  concern,  and  solution  of  the  primary  problem  will 
have  been  completed. 


Table  4.  Kffcct  of  Area  Expansion  Ratio  of  Debris  Pit  on 
Inlet  Duct  I- low  and  Blast  Valve  Duct  Flow1* 


Area 

Kxpansion 

Ratio 

Shock 

Propagation 

Tin'c^ 

(see) 

Parameters  at 
Fxit  to  Inlet  Duef 

Shock  Parameters 
in  Blast  Valve  Duet 

Interface 
Parameters 
in  Blast 
Valve  Duet 

P 

(psia  l 

T 

<"K) 

u 

(cm/sec 
X lo4> 

P 

ipsia) 

T 

(°K> 

u 

(cm  'see 
x HI-4) 

X 

(cm) 

AT 

(”K) 

X 

(cm) 

1.0 

0.001(80 

247 

BSI 

8.27 

82.6 

546 

5.05 

m 

■**31 

317 

5.0 

0.00807 

24*» 

8.44 

81.3 

553 

4.98 

KM 

310 

10.0 

0.00842 

25o 

3.692 

8.44 

81.3 

545 

mm 

295 

16.0 

000881 

234 

3.664 

8.75 

80.9 

543 

Sal 

EB 

ESI 

319 

•'  Configuration  of  Figure  2 with  l-Mt  nuclear  surface  wave  at  1,500-foot  radius. 
h Number  of  computation  cycles  was  I .<HIO  for  all  cases. 


Kulcrian  Computer  Code  Results 

The  l.ulcrian  Code  results  are  given  for  a constant  area  duet  only.  The  variable  area  option  for  this 
code  is  not  accurate  and  should  not  be  used,  l-'or  variable  area  problems  the  l.agrange  Code  is  more  appro- 
priate. Since  the  i.agrange  Code  was  in  agreement  with  experimental  data,  the  output  of  the  Kulcrian  Code 
was  compared  to  the  l.agrange  Code  results  to  verify  the  F.ulerian  Code  accuracy.  This  comparison  is 
presented  in  Figure  1 ..  The  time  waveforms  for  overpressures  are  given  at  four  locations  in  a duct  through 
which  an  attenuating  shock  wave  is  propagating.  The  Hitlerian  solution  agrees  well  with  the  l.agrange  solution 
except  for  an  overshoot  at  the  shock  front.  This  overshoot  characteristic  of  the  l.ax-Wendroff  finite- 
difference  scheme  that  is  used  | 2| . The  shock  front  is  steeper  in  the  Fulerian  solution,  and  the  shock  speed 
is  somewhat  faster  th  in  in  the  l.agrange  solution  This  difference  in  shock  speeil  is  attributed  to  the  pressure 
overshoot  which  affects  .hock  front  attenuation  due  to  rarefaction  decay. 

Studies  concerned  with  reflected  shock  waves  at  a rigid  boundary  show  the  reflected  wave  to  have  the 
characteristic  overshoot.  But  the  reflected  pressure  after  the  overshoot  agrees  with  Kankinc-llugoniot  relation- 
ship predictions. 
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A description  of  the  input  and  output  parameter*  that  are  used  in  the  Eulerian  Code  with  information 
on  obtaining  a Code  listing  is  given  in  Appendix  B. 


CONCLUSIONS  AND  RECOMMENDATIONS 

Both  the  Lagrange  computer  code  and  the  Eulerian  computer  code  presented  arc  shown  to  he 
adequate  for  predicting  shock  wave  propagation  in  a constar, ' area  duct.  For  shock  wave  prediction  in 
a variable  area  duct  (he  Lagrange  Code  should  he  used.  Both  codes  adequately  predict  shock  wave  reflec- 
tion from  a rigid  boundary  and  are  suitable  for  either  constant  strength  shock  waves  or  shock  waves  with 
a rarefaction  region  with  exponential  pressure  decay  behind  the  shock  front.  For  shock  wave  prediction 
in  an  air  entrainment  system  consisting  of  debris  pit  and  blast  valve  ducts  (Figure  2)  the  Lagrange  com- 
puter code  should  be  used,  in  which  case  the  simulated  surface  wave  can  be  of  1 'It  nuclear  explosion 
or  TNT  explosion  origin. 
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Appendix  A 

LAGRANGE  COMPUTER  CODE  DESCRIPTION 


INTRODUCTION 

The  Lagrange  computer  code  is  similar  in  basic  structure  to  the  WUNDY  code  reported  in  Reference  17, 
which  was  basal  on  the  KO-CODE  of  the  University  of  California  Radiation  Laboratory  1 18| . The  basic  struc- 
ture of  the  Lagrange  Code  and  the  input  and  output  formats  are  described  below.  The  computer  code  listing 
and  card  deck  can  be  obtained  by  requesting  a program  tape  from  the  Computer  Center,  Code  1.06,  CEL.  All 
quantities  in  the  code  are  given  in  cgs  units  except  for  printout  of  pressure  whirh  is  given  in  psia. 


CODE  BASIC  STRUCTURE 


The  computer  code  consists  of  a mam  control  subroutine  with  several  auxiliary  subroutines.  A list  of 
these  subroutines  with  a description  of  each  function  follows. 


Subroutine 


Function 


MAIN 

BDY1 

BDY123 

BDY2 

DATEXP 

EQST 

EQS1 

F.QS3 

GENR 

GF.OM 

IITEMP 

IIYDR 

NUBDY 

OUT1 

OUT2 

OUT3 

OUT4 

REZ1 


Controls  main  logical  flow  and  reads  input  data 

Specifies  motion  of  interface  at  duct  inlet  for  low  temperature  wave 

Specifies  motion  of  interfaces  at  T-junction 

Specifics  motion  of  interface  at  a duct  exit 

Data  source  for  duct  inlet  losses 

Controls  equation  of  state  subroutine  acquisition 

Equation  of  state  for  ideal  gas  (T  < t,000°K) 

Equation  of  state  for  real  air  (T  < 24.000°K) 

Initializes  problem 

Calculates  cross-sectional  area  and  zone  volume 
Calculates  Z for  EQS3 
Computes  hydrodynamic  motions 

Specifies  motion  of  interface  at  duct  inlet  for  1-Mt  nuclear  wave  case 
Prints  normal  output,  pressure,  etc.,  in  each  zone  at  fixed  times 
Accumulates  data  on  main  shock  front  in  each  duct 
Accumulates  pressure,  etc.,  versus  time  at  fixed  positions 
Punches  cards  from  which  the  problem  can  be  continued 
Removes  excessively  compressed  zones 
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Subrouting 

KtZENl 

RCZF.Nl 

REZEX1 

viMfc&r 


Function 

Ad  lie  • aono  at  duct  entrance  (or  mast  inflow 

Controls  entrance  tone  site  « T junction 
Controls  exit  zone  arte  at  T-juncrion 
Calculate!  time  step 


CODE  INPUT  QUANTITIES  AND  FORMATS 

The  input  dnu  aymbols,  including  the  data  card  number  on  which  they  appear,  and  the  dm  format 
arc  given  in  Table  a-L  The  data  appear s on  a particular  card  in  the  order  in  which  it  it  given.  A sample  data 
card  Hiring  i.l  given  in  Table  A-2  to r a l.JtHHoot  radius  location  from  a t-Mt  nuclear  burst  with  an  ait  entum- 
nieot  geometry  as  shown  in  Figure  12. 

CODE  OUTPUT  VARIABLES  AND  FORMATS 

The  output  of  this  computer  code  consists  of  printout  of  the  input  dm.  corrected  input  data.  ouTl 
subroutine  printout.  OLT2  subroutine  printout,  and  0UT3  mfaioutinc  printout.  The  OUT*  subroutine  punches 
cards.  To  give  a full  output  would  be  too  lengthy:  therefore,  only  a sample  output  from  printout  of  rh<  input 
data  and  Ot/Tl  at  cycle  500  are  given  5n  T«bl«*  A~l  and  A-4-,  r«*p*etiv«1y.  The  data  in  Table*  A.2  rh rough  A>4 
are  for  the  problem  of  Figure  12.  The  program  was  run  On  a CDC  0600  computer,  and  it  reqvired  a core  storage 
of  105.000  • M a running  rime  of  1 19  seconds. 

The  normal  output  is  provided  bp  the  0UT1  subroutine,  which  prints  out  the  velocity,  displacement, 
and  several  state  variables  for  each  zone  at  desired  times,  The  printout  is  controlled  by  the  quantity  NPR  me 
variable*  printed  our  every  NPR  cycles  of  computation  arc  given  in  Table  A-5. 

An  auxiliary  output  is  provided  by  the  OUTJ  sufaroviioa,  which  print*  out  variables  at  desired  locations 
in  the  duct  versus  time.  The  control  variable  5(1, J)  specifics  the  location  at  which  the  variable*  given  in  'fable  t\  6 
are  printed  out. 

Table  A-l.  Jnpur  Quantities 

Catd 


Number 

Format 

Symbol 

Definition 

1 

7A10 

ALiST(l) 

Problem  identification 

2 

7A10 

ALIST(l) 

Problem  identification  continued 

3 

1215 

NPROB 

Problem  number 

IMAXL 

Tocal  number  of  ducts.  1 or  3 

intape 

* 0,  no  data  input  from  tape  18 

tNCOPS 

* O,  no  extra  input  from  card* 

continued 
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Table  A-l.  Continued 


Card 

Number 


4 


5 


6 


7 


Format 

Symbol 

Definition 

NQU1T 

Total  number  of  cycles  to  run 

NPR 

Print  after  every  NPR  cycle 

NTAPE 

■ 0,  do  not  write  tape  1 8 

KOPT 

» 1,  side-on  type  entrance 
■ 0.  wave  originates  at  duct  entrance 

1213 

KOUT2 

■ 0,  do  not  call  OUT2  (usually  4 when 
OUT2  is  used) 

KOUT2A 

Store  data  every  KOUT2  eydes 

KOUT2B 

Controls  coupling  between  OUT2  and 
OUT3i  usually  zero 

KOUT3A 

Store  P-T  data  every  KOUT3  cycles 

KOUT4 

Punch  continuation  cards  at  KOUT4 
cycles 

KRKZ1 

* 0.  REZI  not  used 

- 1,  REZI  used 

1215 

NC(1) 

through 

NC(12) 

Control  variables,  all  zero  except  NC(6> 
and  NC(7) 

NC(6) 

* 1 , print  NC(6)  times  per  decade  in  time 

NC(7> 

* 1 . use  special  pseudw-viscoiity 

1215 

NC(13> 

through 

NC(24) 

Control  variables,  all  zero  except  NC(  17) 

NC(  1 7) 

» 1.  print  CRAMS  instead  of  PDVN 

7E10.4 

EBI 

* 0,  not  now  used 

T 

= 0,  start  with  time  zero 

DTMIN(2) 

“ 0,  use  built-in  time  step 

DTRATE 

= 0,  use  built-in  time  step  change  rate  of 
1.4 

STABIL 

= 0,  use  built-in  stability  constant  of  0.81 

UCUT 

= 0,  use  built-in  velocity  cutoff  value  of 
1 x ltr2  cm/sec 

continued 
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Table  A- 1 . Continued 


Card 


Number 

bormat 

Symbol 

Definition 

Surface  wave  constants 

H 

r;:io.4 

Al 

through 

A5 

U 

7i:io.4 

HI 

through 

K7 

Surface  wave  exponents 

10 

7 1: 10.4 

1)1* 

Positive  pressure  phase  duration 

DU 

Positive  velocity  phase  duration 

IS 

Shock  arrival  time 

TH 

Time  of  maximum  temperature 

PSO 

Initial  value,  nuclear  wave  pressure 

(JSO 

Initial  value,  nuclear  wave  dynamic 
pressure 

II 

7KKI.4 

IMPS  1 

Surface  temperature  constant 

T.MPS2 

Surface  temperature  constant 

II.AC 

s 1 , use  NUBDY  subroutine 

“ use  BDY1  subroutine 

12 

7i-:io.4 

TI.ISTd) 

through 

ri.lST(A) 

All  zero,  printing  controlled  by  NPR 

13 

515 

1 = 1 

Duct  identification 

NKQSTd) 

Number  of  equation  of  state  in  duct  1 

jcAi.cn  > 

Number  of  last  interface  currently  being 
calculated  in  duct  1 

NZONKSd) 

lotal  number  of  zones  in  duct  1 

KOUTJd) 

Store  P- 1 data  at  KOUT3  locations  in 
duct  1 

14 

HK  If  Ml 

CAMMAId) 

(iamma  used  in  duct  1 

OUTBDYdl 

l.ength  of  duct  1 

l-INIKI) 

Initial  internal  energy  in  duct  1 

UINlTd) 

Initial  velocity  in  duct  1 

DINIKl) 

Initial  density  in  duct  1 

continued 
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Table  A- 1.  Continued 


Card 


Number 

Format 

Symbol 

Definition 

PRICTd) 

Friction  factor  in  duct  1 

CINQd) 

2.0,  pseudo-viscosity  constant,  duct  1 

AINQd) 

0.2,  pseudo-viscosity  constant,  duct  1 

15 

6E10.0 

DOd) 

Diameter  for  Xd.JXXId) 

Dl(l) 

Linear  rate  of  change  in  diameter 

D2(l) 

■ 0,  not  used 

Did) 

■ 0,  not  used 

xid) 

Begin  linear  diameter  change  at  Xid ), 
duct  1 

X2(l ) 

End  linear  diameter  change  at  X2(  1 ),  duct  1 

16 

6E10.0 

Sd.l) 

through 

S(l,6) 

Positions  in  duct  1 to  rollcct  P-T  data 
by  OUT3;  zero  if  not  used 

!:  Repeat  cards  15  through  16  for  1 

= 2.IMAXL. 

Last  card 

IS 

NEXT 

* 1,  read  new  set  of  data 

=♦  I,  stop;  end  of  computation 


Table  A-2.  Data  Card  Listing- Problem  of  Figure  12. 
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U«* 
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3. 
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0.0 

-r  1 

/.<* 

• 00 
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♦ 02 

1.5 

•02 

3.S 

♦ 02 
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♦ 00 
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♦ 00 
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♦ 00 
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1. 

• 1A 

r. 

♦00 

0. 

• or 
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0. 
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Table  A* 3.  Printout  of  Input  Data  Problem  of  Figure  12. 
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Tabic  A-4.  Continued 
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Table  A-J.  OUT!  Subroutine  Printout  Variable* 


Printout  Variable  De*cription 


IDT 

Duct  with  tmallest  time  step 

JDT 

Zone  with  smallest  time  step 

PS 

Surface  pressure 

QS 

Surface  dynamic  pressure 

TEMPS 

Surface  temperature 

ZS 

Surface  compressibility  factor 

GAMMAS 

Surface  value  of  gamma 

HTS 

Surface  value  of  total  enthralpy 

DS 

Surface  density 

ES 

Surface  internal  energy 

1 

Number  of  duct 

J 

Number  of  interface 

XO.J> 

Distance  to  interface 

U(I.J) 

Velocity  of  interface 

XAV(I.J) 

Distance  to  zone  center 

DO.J) 

Density 

E(I.J) 

Internal  energy 

PQO.J) 

PtQ 

QG.J) 

Pseudo-viscosity 

PU.J) 

Pressure 

GRAMS(I.J) 

Zone  mass 

TEMPO,  J> 

Temperature 

DIA(I.J) 

Duct  diameter 

OTZJ(IJ) 

Time  step  in  duct  1 at  zone  ) 

GAMMAO.J) 

Adiabatic  exponent  gamma 

ZJ<I,J) 

Compressibility  factor 

Table  A-6.  OUTJ  Subroutine  Printout  Variable* 


Printout  Variable* 

Description 

1 

Duct  number 

S(I.J) 

Position*  P-T  data  collected 

ncyci.k 

Number  of  computation  cycle* 

T 

Real  problem  time 

PSI(I.J) 

Pressure,  psia 

OVPSI(I.J) 

Overpressure,  psi 

PDYN(I.J) 

Dynamic  pressure,  psi 

rxi.j) 

Density 

Ud.J) 

Velocity  of  interface 

TEMPO, J) 

Temperature 
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Appendix  B 

EULERIAN  COMPUTER  CODL  INPUT  AND  OUTPUT  PARAMETERS 


The  input  anti  output  parameters  are  presented  for  the  Eulerian  computer  code.  A tape  of  the  code 
for  listing  and  card  punching  is  available  upon  request  from  the  Computer  Center,  Code  1.06,  CEL. 


COMPUTER  CODE  INPUT  PARAMETERS 

The  input  parameters  are  listed  in  Table  B-l. 


Table  B-l.  Eulerian  Code  Input  Parameters 


Symbol 

Definition 

PO 

Ambient  pressure,  lb/in.2 

TO 

Ambient  temperature,  °R 

RO 

Ambient  specific  gas  constant,  fr/°F 

CO 

Ambient  specific  heat  ratio 

R2 

High  temperature  specific  gas  constant,  ft/°F 

(12 

High  temperature  specific  heat  ratio 

Al(i),  A2(i),  A3(i), 
A-Ki):  i - 1.4 

Constants  in  equation  for  duct  radius  in  duct  section  i, 
r(i)  * Al(i)  + A2(i)x  + A3(i)x2  + A4<i)  x3 : all  zero 
except  Al(l)  for  fisting  given  in  Appendix  B. 

XB(  1 ),  XB(2),  XB< 3 ) 

Maximum  x values  in  the  1st,  2nd,  and  3rd  sections  of 
the  duct;  all  greater  than  XFIN  for  listing  given  in 
Appendix  B. 

Bl,  B2,  B3,  B4 

Constants  in  the  equation  for  overpressure  transmitted 
into  the  duct;  p<0,0)  ■ pa  + Bl(pa  - pa)B2;  (B3  and  B4 
arc  not  currently  used) 

XIN 

X value  at  duct  entrance;  use  XIN  * 0 

XHN 

X value  at  end  of  duct,  ft 

DEI.X 

Ax,  finite-difference  mesh  size,  ft 

A 

ALPHA 

Mach  number  of  nuclear  blasr  wave  as  it  passes  over 
duct  entrance 

— t/t ; 

tj,  time  constant;  p(0,t)  = pa  + |p(0,0)-pa|e 

BETA 

Not  used 

SI.S2 

Constants  used  in  computing  friction  factor,  f; 
f = Sl(Re)s2 

continued 
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Table  B-2.  Continued 


Symbol 


Definition 


XXI Y 
Kl,  K2,  K3 
ID 

IY 

NEND 

NTIME 

IWRTTP 

IWRTXP 

NF 

in. 


NWRT 

IWP 


K,  constant  used  in  computing  time  increment; 

At  « Ax/F(U  + C) 

Constants  used  in  the  computation  of  the  pseudo* 
viscosity  factors,  usually  2.0 

• 1,  if  duct  is  rectangle  of  unit  width 

■ 2,  if  duct  is  circular 

Not  used 

Initializing  parameter,  set  NKND  * I 

Maximum  value  of  n,  i.e„  number  of  time  steps  to  be 
computed 

Output  is  printed  every  IWRTTP,h  time  step 

Output  is  printed  every  IWPTXP**1  space  step  for  every 
IWRT1'P,h  time  step 

» 0,  if  friction  factor,  f,  is  zero 

c 2 

« I,  if  friction  factor,  f,  is  computed  from  f » Sl(Re) 

Number  of  time  steps  that  must  l>e  computed  before  the 
pseudo-viscosity  factors  take  on  their  full  value;  for 
n < IT!.,  the  pseudo-viscosity  factors  are  proportional 
to  n/ITI. 

The  maximum  value  of  pressure  at  a given  value  of  axial 
distance;  prinred  out  every  NWRT,h  space  step 

Properties  behind  the  nuclear  blast  wave;  printed  out 
every  iWPlh  time  step 


COMPUTER  CODE  OUTPUT  PARAMETERS 

The  output  parameters  are  listed  in  Table  B-2. 


Table  B-2.  Kulerian  Code  Output  Parameters 

Symbol  Definition 

TIME  Elapsed  time  from  appearance  of  shock  wave  at  duct 

entrance,  sec 

X Axial  distance,  ft 

continued 
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Table  R-2.  Continued 


Symbol 


Definition 


X/l. 

KUO 

Y 

Y/Yl 

Li 

AKKA 

ARKA/ARKAl 

MACH  NO 
P 

P/PO 

RIIO/RIIO 

T 

T/TO 

I 

X(l) 

X/D 

P/PO 

P 


Nondimcnsional  axial  distance,  normalized  with  respect 
to  duct  length 

Density,  lbm/ft5 

Duct  radius  of  half-height,  ft 

Nondimcnsional  duct  radius  or  half-height,  normalized 
with  respect  to  initial  value 

Velocity,  ft/sec 

Duct  area,  ft2 

Nondimcnsional  duct  area,  normalized  with  respect  to 
initial  value 

Mach  number 

Pressure,  lb/in.2 

Nondimcnsional  pressure,  normalized  with  respect  to 
ambient  value 

Nondimcnsional  density,  normalized  with  respect  to 
ambient  value 

Temperature,  °F 

Nondimcnsional  temperature,  normalized  with  respect  to 
ambient  value 

Number  of  space  increment 

Value  of  X at  the  Ith  space  increment,  ft 

Nondimcnsional  axial  distance,  normalized  with  respect  to 
initial  duct  diameter  or  height 

Maximum  value  of  nondimcnsional  pressure  at  X(l), 
normalized  with  respect  to  ambient  value 

Maximum  value  of  pressure  at  X(l),  lb/in.2 


BASIC  NOMENCLATURE 


A 

A,-As 

a 


b,-b5 


D 


o 


♦ 

u 


Duct  cross-sectional  area 
Constants,  Lagrange  Code 
Acceleration 


Psen  lo-viscosity  constants,  Lagrange  Code 
Constants,  Lagrange  Code 
Acoustic  velocity 


Duct  diameter 

Duration  of  positive  pressure  phase,  nuclear 
wave 


Duration  of  positive  velocity  phase,  nuclear 
wave 


E 

F 

F<U) 

f 

»2 

*T 

G(U) 

h 


Specific  internal  energy 
Variable  defined  aspte  ♦ l/2u2) 
Stability  constant.  Eulcrian  Code 
Function  of  vector  U 
Friction  factor 

Pseudo-viscosity  factor,  F:ulerian  Code 
Pseudo-viscosity  factor,  Eulerian  Code 
Pscudc  viscosity  factor,  Kulerian  Code 
Pseudo-viscosity  factor,  F.ulcrian  Code 
Pseudo-viscosity  factor,  Eulcrian  Code 
Pscudo-viscositv  factor.  Eulerian  Code 
Function  of  vector  U 
Specific  enthalpy 
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Total  enthalpy  per  unit  matt 
Duct  hydraulic  radius 
Pseudo-viscosity  constant,  Kulerian  Code 
Duct  length 

Lagrange  variable,  mass  per  zone 
Mass  change  per  zone.  Lagrange  Code 
Shock  wave  Mach  number 
Variable  defined  as  p u A,  Kulerian  Code 
Variable  defined  as  p A,  Kulerian  Code 
Pressure 

Pressure  at  position  x at  time  t 
Variable  defined  as  p ♦ q,  Lagrange  Code 
Pseudo-visa  -,itv,  Lagrange  Code 
Dynamic  pressure 
Particular  gas  constant 
Reynolds  numtier 

Friction  factor  constants,  Kulerian  Code 

Specific  entropy 

Time 

Shock  arrival  time 
Initial  slope  time  intercept 
Time  increment 
Particle  velocity 
Shock  wave  velocity 


Variable  vector,  Kulerian  Code 


V 

x 

Ax 

Z 


*,.02.0, 


7 

f 


P 

P<x,t> 


r 


T 


w 


U) 


Volume  per  unit  nun 

Distance  along  duct,  Kulerian  variable 

Distance  increment 

Compressibility  factor  for  real  gas 

Pseudo-viscosity  parameter,  Eulerian  Code 

Adiabatic  exponent 

Small  quantity 

Defined  as  <7  - I )/(y  + , ) 

Density 

Density  at  position  x at  time  t 
Defined  as(t-t,)/D* 

Shear  stress  at  wall  of  duct 
Defined  as  (t  - t,)/D* 


Subscripts 

* Downstream  of  shock 

2 Upstream  of  shock 

* Ambient  condition  value 

e Duct  inlet  value 

J Net  point  location 

0 Initial  value 

* Surface  value 


Superscripts 


n 


Number  of  time  increments 


